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ABSTRACT 

Depending on mass and metallicity as well as evolutionary phase, stars occasionally experience convective- 
reactive nucleosynthesis episodes. We specifically investigate the situation when nucleosynthetically unprocessed, 
H-rich material is convectively mixed with a He-burning zone, for example in convectively unstable shell on top 
of electron-degenerate cores in AGB stars, young white dwarfs or X-ray bursting neutron stars. Such episodes 
are frequently encountered in stellar evolution models of stars of extremely low or zero metal content, such as 
the first stars. We have carried out detailed nucleosynthesis simulations based on stellar evolution models and 
informed by hydrodynamic simulations. We focus on the convective-reactive episode in the very-late thermal 



pulse star Sakurai's object (V4334 Sagittarii). |Asplund et al. ( 1999) determined the abundances of 28 elements 



many of which are highly non-solar, ranging from H, He and Li all the way to Ba and La, plus the C isotopic ratio. 
Our simulations show that the mixing evolution according to standard, one-dimensional stellar evolution models 
implies neutron densities in the He intershell (< few 10 n cm" 3 ) that are too low to obtain a significant neutron 
capture nucleosynthesis on the heavy elements. We have carried out 3D hydrodynamic He-shell flash convection 
simulations in 4ir geometry to study the entrainment of H-rich material. Guided by these simulations we assume 
that the ingestion process of H into the He-shell convection zone leads only after some delay time to a sufficient 
entropy barrier that splits the convection zone into the original one driven by He-burning and a new one driven 
by the rapid burning of ingested H. By making such mixing assumptions that are motivated by our hydrodynamic 
simulations we obtain significantly higher neutron densities (~ few 10 15 cm" 3 ) and reproduce the key observed 
abundance trends found in Sakurai's object. These include an overproduction of Rb, Sr and Y by about 2 orders of 
magnitude higher than the overproduction of Ba and La. Such a peculiar nucleosynthesis signature is impossible to 
obtain with the mixing predictions in our one-dimensional stellar evolution models. The simulated Li abundance 
and the isotopic ratio 12 C/ 13 C are as well in agreement with observations. Details of the observed heavy element 
abundances can be used as a sensitive diagnostic tool for the neutron density, for the neutron exposure and, in 
general, for the physics of the convective-reactive phases in stellar evolution. For example, the high elemental 
ratio Sc/Ca and the high Sc production indicate high neutron densities. The diagnostic value of such abundance 
markers depends on uncertain nuclear physics input. We determine how our results depend on uncertainties of 
nuclear reaction rates, for example for the 13 C(a,n) 16 reaction. 

Subject headings: stars: AGB and post-AGB — stars: abundances — stars: evolution — stars: interior — stars: 
individual (V4334 Sagittarii) — physical data and processes: hydrodynamics — physical data 
and processes: nuclear reactions, nucleosynthesis, abundances 

1. INTRODUCTION He-shell flash convection of Asymp totic Giant Branch ( AGB) 

, , „ .. , r . „ , .. stars, such as the branching at 128 I (Reifarth et al.|2004i. This 
1.1. Convective-reactive phases of stellar evolution ... . , °. > — ■. , ■ , 



situation is comparatively simple to simulate as the rapid nu- 
In stellar evolution the nuclear time scale is usually much clear reaction in question, the double-decay of m l, does not 
larger than the convective mixing time scale. However, this release any significant amount of energy. A post-processing ap- 
is not always the case. An example of stellar nucleosynthe- proach of the standard stellar evolution calculation with some 
sis where nuclear reactions and convective mixing occurs on one-dimensional treatment of convection, like mixing-length 
the same time scale are slow neutron capture process branc h- theory (MLT), with time-dependent mixing gives a reasonable 



ings (s process, Burbidge et al.|1957||Wallerstein et al.|1997] > in approximation of this situation^] 
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10 Although even in this case multi-dimensional effects of convection have to be taken into account eventually as simulations by Herwig et al 
velocity profile at the bottom of the convective shell is flatter compared to the MLT prediction. 
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The goal of this paper is instead to investigate the situation 
when rapid nuclear reactions are indeed releasing amounts of 
energy that are likely to affect the fluid flow, as for example in 
the case of proton capture of 12 C in convective He-burning. In 
the fluid dynamics community this mixing regime is sometimes 
refered to as level-3 mixing, where the flow is altered by atten- 
dant changes in the fluid (Dimotakis 2005). We refer to these 
situations as reactive-convective phases in order to emphasize 
the fact that the time scales of highly exothermic nuclear reac- 
tion and the convective fluid flow time scales are of the same 
order. 

The ratio of the mixing time scale and the reaction time scale 
is called the Damkohler number: 

D Q = — . (1) 

7"react 

MLT is concerned with averaged properties both in time over 
many convective turn-overs and in space over the order of a 
pressure scale height. In the categories of |Dimotakis| (|2005 ) 
diffusion coeffiecients derived from MLT may describe level- 
1 mixing (while mixing induced by rotation involves flow dy- 
namics that are altered by mixing processes and labeled in this 
scheme as level-2 mixing). Therefore, time-dependent mixing 
through a diffusion algorithm with diffision coefficients derived 
from MLT is appropriate for regimes with D a <C 1. The diffi- 
culty of simulating convective-reactive phases in present one- 
dimensional stellar evolution codes then appears as the inability 
of MLT (or any similar convection theory) to properly account 
for the additional dynamic effects introduced through rapid and 
dynamically relevant nuclear energy release in level-3 mixing 
associated with Damkohler numbers D a « 1 . 

Convective-reactive episodes can be encountered in numer- 
ous phases of stellar evolution, including the He-shell flash 



of AGB stars of extremely low metal content (e.g. Fujimoto 


et al. 


2000 


|Suda et al.||2004| |Iwamoto et al.||2004| Cristallo 


et al. 


2009 


i, metallicity low-mass stars (e.g. [Hollowell et al. 


1990||Schlattl et al.|2002||Campbell & Lattanzio 
white dwarfs of solar metallicity (e.g. Iben et ak] 


|2008|), young 
1983||Herwig 



rotating Pop III massive stars (Ekstro m et al.|2008]|, and more 



in general, in low metallicity massive stars (Woosley & Weaver 
1995). These combustion events are encountered as well in 



X-ray burst calculations of accreting neutron stars (Woosley 
|et al.|2004||Piro & Bildsten|2007| ), and accreting white dwarfs 



(Cassisi et al. 1998 1 that may be the progenitors of SN la. 
Convective-reactive events have been found in post-RGB stel- 
lar evolution models and associated with the horizontal branch 



anomalies in certain globular clusters ( Brown et al. 2001 Miller 
Bertolami et al. 2008). Finally, again in AGB stars, convective- 
reactive phases can be found in hot dredge-up ([Herwig 2004; 
|Goriely & Siess|2004| |Woodward etaTpOQg) , a phenomenon 
that is associated with the treatment of convective boundaries, 
generally in more massive and lower metallicity AGB stars. 

Although convective-reactive phases are quite common in 
stellar evolution, in particular in the early, low-metellicity Uni- 
verse, we do not currently have a reliable and accurate way of 
simulating them. In this work we discuss the case of the He- 
shell flash with H-ingestion in a very-late (post-AGB) thermal 
pulse at solar metallicity. This situation is extremely similar to 
H-ingestion associated with the He-shell flash in AGB stars at 
extremely low metallicity. The one-dimensional, spherically- 
symmetric stellar evolution approximation is not very realis- 
tic in this case, because both the entrainment of H into the 
He-shell flash convection zone as well as the subsequent con- 



vective transport, mixing and nuclear burning of hydrogen en- 
riched fluid parcels are inherently a three-dimensional hydro- 
dynamic process. The energy from proton captures by 12 C via 
the 12 C(p,7) 13 N reactions is released on the same time scale 
(~ 1 . . . lOmin for T = 1.3 ... 1 .05 x 10 8 K) of the fluid flow of 
convection (Sect.|B|, and this energy will add entropy to fluid 



elements and in turn feedback into the hydrodynamics (Her- 
wig 2001 ). These highly coupled, multi-dimensional processes 
are approximated in through the MLT, complemented with a 
time-dependent mixing algorithm. This assumption may not be 
realistic in the present case (see Sect. 3.2 and 4.2 1. 



1.2. Post-AGB flash star Sakurai's object and its observed 
abundance properties 

Sakurai's object is a very-late thermal pulse post-AGB ob- 
ject (Duerbec k et al.|2000 and ref. there) and has experienced 
a H-ingestion flash in 1994. The star's observed abundance sig- 
natures are highly non-solar, and very unusual for a post-AGB 
low-mass star (Sect. |3.2j ). Nevertheless, there is wide agreement 
in the literature that the object's distance is 2 - 5kpc and that it 
has a mass of around 0.6 M (van Hoof et al. |2007| and ref. 
therein), pointing to a low mass star progenitor. Moreover, the 
high abundance of Li requires the existence of 3 He in the enve- 
lope (Herwig & Langer 2001), pointing again to a low mass star 
progenitor that was not affected by hot bottom burning (HBB). 
Indeed, hot-bottom burning occurs at solar metallicity for stars 
with Mzams ~ 4M Q and destroys 3 He in the AGB envelope 
(Scal o et al.|1975| l. Another process, that could effect the evo- 
lution of J He during the progenitor evolution of Sakurai's object 
is extra-mixing below the convective envelope during either the 
RGB or AGB (e.g. |Wasserburg et al.|1995||Charbonnel & Zahn| 
2007 1 Denissenkov||2010| i. Sakurai's object serves potentially 



as an important constraint for theories of such mixing because 
the observed Li abundance increase during the observations in 
1996 as reported by Asplund et al.| ( |T999] > can only be modeled 
in the very late thermal pulse if significant amounts of 3 He are 
still present in the envelope at the beginning of the post-AGB 
evolution. 

The light curve of this object was closely monitored as it 
evolved within approximately 2yr from the pre- WD location 
in the HRD back to the AGB location, a much shorter evolu- 
tion time scale than previously predicted ( |Herwig et al.|l9 99). 
A possible explanation of such a fast born-again evolution of 
Sakurai's object is that the convective mixing efficiency in the 
He-shell flash convection zone is smaller by a factor of ~ 30 
compared to the MLT predictions in standard one-dimensional 
stellar models (Herwig 2001). This modification is motivated 
by the reasoning that in the convective-reactive regime the fluid 
flow would be eventually strongly affected by the energy re- 
leased rapidly on a time scale comparable to the fluid flow ve- 
locity. This process, indeed, would locally add boyuancy to 
the fluid element causing a behavior that is not reflected in the 
mixing-length theory. 

Miller Bertola mi et al.| ( |2006| l have presented a more de- 
tailed investigation and emphasize the importance of appropri- 
ate time resolution. In addition, they studied the role of over- 
shooting and /i-gradients. Their simulations with exponential, 
depth-dependent overshooting agree better with observations 
than tracks computed without any overshooting, /i-gradients 
appear to have only secondary effects. Confirming the mass 
dependence of the proton-ingestion born-again evolution first 
reported by |Herwig| ( |2001| l, |Miller Bertolami & Althaus| ( |2007| l 
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point out that the initial return light-curve of Sakurai's object 
could be reproduced with a slightly lower mass model than the 
0.604M Q adopted by Herwig (2001), a high time resolution 
and their alternative description of convective transport. How- 
ever, the second heating phase into which the Sakurai's object 
has entered now ( |van Hoof et al.||2007| ), seems to be better in 
agreement with the modified convection models proposed by 
|Herwig| l |200T) . 

While the light curve of Sakurai's object has certainly raised 
doubts about the capability of one-dimensional stellar evolu- 
tion calculations to reproduce its evolution, in this work we 



show that the abundance determinations by Asplund et al. 



(1999] pose a much more stringent constraint on the physics 
of convective-reactive phases. Asplund et al. determined 28 el- 
emental abundances at four times between April and October 
1996, when the star had cooled to below 8000 K. In particu- 
lar, among light elements a significant enhancement (at least 
0.5 dex) with respect to the solar abundance has been observed 
for Li, Ne and P. Beyond iron, Cu, Zn, Rb and Sr peak ele- 
ments are significantly enhanced. In addition, there are trends 
as a function of time that are smaller than the differences to so- 
lar. However, for this initial analysis which is not yet based on 
full hydrodynamic simulations with nuclear burn, we will not 
discuss those trends in detail. 

A few preliminary comments on individual elements may 
be in order. The observed Li is clearly produced above the 
meteoritic value. |Herwig & Langer] ( |2001| l proposed that to- 
gether with protons J He is ingested into the He-shell flash con- 
vection zone, providing the fuel to produce Li via the reac- 
tion chain 3 He(a,7) 7 Be(/3 + ) 7 Li. The first s-process peak ele- 
ments are enhanced by up to 2dex while Ba and La are not 
enhanced, causing a ratio of Ba peak to Sr peak elements 
that is much lower than expected from models and observa- 
tions of AGB stars ( |Busso et aT]|2001[) 



We can translate the 



abundances observed by Asplund et al.| ( | 1999[ > into the ra- 
tio of the two s-process indicator indices hs and Is. An s- 
process index s/sq is the overproduction factor of a group of 
i-process elements with respect to the initial solar value. The 
index ratio [hs/ls] = [hs/Fe] - [ls/Fe] monitors the distribution 
of the ^-process elements, and it is an intrinsic index of the 
neutron ca pture nucleosynthesis on heavy elements (Luck & 



Bond 



19911. We have used [ls/Fe]: 



4([Sr/Fe]+[Y/Fe]+[Zr/Fe]) 



and [hs/Fe]=j([Ba/Fe]+[La/Fe]), where square brackets indi- 
cate the logarithmic ratio with respect to the solar ratio (Ta- 
ble[T|. For Asplund's October measurements the indices are 
[hs/Fe]= 0.05 and [ls/Fe]= 1.9 assuming that [Fe/H]= 0.0 for 
Sakurai's object. We record measurements of ±0.2 . . .0.3dex 
as the average approximate index ratio [hs/ls] ~ -2 at the end 
of the observed period. In Fig.[T| we compare such ratio with 
i-process theoretical predictions and stellar observations of low 
mass AGB stars, that are the progenitor population of the Saku- 
rai's object. In particular, we show that the observed [hs/ls] is 
a factor of ten or more lower than in typical AGB stars. There- 
fore, the nucleosynthesis environment that has generated the 
abundances observed by Asplund et al. was very different from 
that encountered in the previous AGB phase. In Fig. [T] we also 
include [hs/ls] from our nucleosynthesis calculations presented 
in this paper, that succesfully reproduce the same ratio mea- 
sured in the Sakurai's object. Such calculations will be dis- 
cussed in detail in Section|5] 

The abundance pattern of Sakurai's object further distin- 
guishes itself from the AGB stars through the significantly en- 



hanced P, Cu and Zn. These elements are not usually produced 
in low-mass stars. Several other elements are reduced, i.e., S, 
Ti, Cr and Fe. In particular, Fe is expected to be depleted, since 
it is the seed for n-capture nucleosynthesis. All these abun- 
dance signatures appear to be the result of a n-capture burst of 
large n-density. Another important feature is the C isotopic ra- 
tio 12 C/ I3 C ~ 3-4, where the large 13 C abundance results from 
the 12 C(p,7) 13 N(/3 + ) 13 C reaction channel. 13 C is also the main 
neutron source during the H ingestion event, which causes the 
peculiar abundance signature observed by Asplund et al. (see 
Section[5]for details). 

In the following we will briefly describe the tools we use 
in this investigation (Sect. [2]) and defer more details to an ap- 
pendix (Sect. [A}. Next we describe the stellar evolution picture 
of Sakurai's object and show how nucleosynthesis simulations 
based directly on the output of one-dimensional stellar evolu- 
tion calculations fail to account for the observed abundance 
patterns (Sect.|3]l. Then we describe hydrodynamic simula- 
tions of entrainment into He-shell flash convection that moti- 
vate our modified mixing assumptions (Sect.|4)>. We will show 
how corresponding nucleosynthesis simulations account for the 
observed abundances, and we discuss the incfluence of nuclear 
reaction rate uncertainty (Sect. [5}. The paper ends with a sum- 
mary and some remarks on implications for the nucleosynthe- 
sis in the first generations of stars, including the light-element 
primary process (Sect.|6]l. In the appendix we give additional 
information on the codes we have used, and in Appendix [B] we 
discuss time scales for burning and mixing. 

2. SIMULATION CODES 

Three different types of simulation codes have been used in 
this work: 

• a stellar evolution code (EVOL), providing one- 
dimensional stellar evolution up to the post-AGBn and 
thermodynamic structures for the beginning of the post- 
AGB He-shell flash event, also known as the very late 
thermal pulse (VLTP); 

• a multi-zone post-processing nucleosynthesis code 
(PPN) with complete nuclear network and mixing; 

• a multi-dimensional-hydrodynamical code (PPM), to 
study how hydrogen is ingested during the VLTP. 

We have used the stellar evolution code EVOL to calculate 

experiencing 



the global evolution of post- AGB stars (Sect.|3.1 
a VLTP ( |BTo^er|1995|perwig|2000[[Herwig & Austin|2004| l 



The assumptions and input physics choices are very similara to 
those in Herwig (2001). Furthermore, we have used structures 
from the last thermal pulse of the AGB model by|Herwig &| 



Austin ( 2004 1, and of the VLTP model by Herwig et al. 



For the detailed nucleosynthesis simulations (Sect. 3.2 



1999). 
_and|5]l 

we have used the PPN (Post-Processing Nucleosynthesis) code 
( |Herwig et al.||2008] >. This code allows to calculate the com- 
plete nucleosynthesis along the radial profile of a star according 
to the structure input from a stellar evolution model in as many 
zones as required. Nuclear burn steps are alternated with time- 
dependent mixing steps. Details, including the nuclear physics 
data information, are given in Sect. 



A.l 



In order to investigate the hydrodynamic behaviour of unpro- 
cessed H-rich material entrained into the He-shell flash convec- 
tion (Sect.Q, we used Woodward's PPM gas dynamics code 
with the PPB advection scheme on a cartesian grid (Woodward 
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et al.|2003[|2006l [Wo odwar d et al.|2008) . For important code 
details, see Sect. |A.2[ 

3. THE STELLAR EVOLUTION PICTURE 

3.1. Global stellar evolution scenario and calculation 

The VLTP evolution scenario involves a He-shell flash on 
a single young white dwarf after the end of H-shell burning 
when the evolution track has just entered the white dwarf cool- 
ing curve in the HRD, as for example shown in Herwig et al. 
(1999}, and in more detail in Sect. 3.2.1 of Miller Bertolami 
et al. ( 2006| l. It involves the convective ingestion of all or 
parts of the small (~ 1O~ 4 M ) remaining unprocessed, and 
thus H-rich, envelope into the hot (T = 1 . . .3 x 10 8 K) He- 
burning flash layers. This He-burning convection zone contains 
a mass fraction of 20 . . . 40% (depending on convective model 
assumptions, Herwig 2000 ; Miller Bertolami et al. 2006 ) of pri- 
mary 12 C. Protons are rapidly captured by the abundant I2 C, 
on the time scale of convective fluid flows of approximately 
5 ... lOmin. 

The progenitor is a low mass AGB star for which s-process 
element enhancements are expected at the Sr-Y-Zr peak and at 
the Ba peak (e.g. Busso et al. 2001 1. elements signature ob- 



served in Sakurai's object is not typical of the s process in AGB 
stars. Indeed, according to the observations by Asplund et al. 
( |1999| l, the ratio of the second peak to the first peak s-process 
elements is [Ba/Y] ~ -2, in contrast to the expected AGB stars 
ratio -1 < [hs/ls] < 1 at solar-like metallicity (e.g. |Busso et al. 
2001). This result does not change if we assume a lower than 
solar metallicity for Sakurai's object of [Fe/H] = -0.63 (values 
between brackets in Tab[T]i. Such a choice may be indicated by 
the sub-solar observed Ba abundance, and indeed, the Ba and 
La abundance even lead us to assume that there was no signifi- 
cant i-process contribution in the previous AGB phase at all. 

In any case, the peculiar abundance signatures of Sakurai's 
object has to originate in the H-ingestion event of the VLTP, 
and can not be explained in terms of any nucleosynthesis dur- 
ing the AGB progenitor evolution. 

The initial abundance distribution for our post- AGB He-shell 
flash nucleosynthesis simulations is a combination of light ele- 
ments (with A < 23) from the intershell abundance of an AGB 
star at the end of the evolution taken from a 2M simulation 
similar to those in Herwig & Austin (2004), and heavier species 
according to Asplund et al. (2005 ) with the isotopic ratios from 
Lodders|([2003|l scaled tometallicity [Fe/H] = -0.18. 



The intershell abundances that matter for our simulations are 
mostly primary He-burning products, so details of the initial 
abundance are not important. The choice of more recent solar 
abundances (Aspl und et al.|[2009| |Lodders et al.||2009) > would 
not modify the results presented in this paper. 

In the following section we will discuss the nucleosynthesis 
according to one-dimensional stellar evolution mixing predi- 
tions of the very-late thermal pulse. 

3.2. Nucleosynthesis according to the stellar evolution model 

Fig. [2] shows the H-profile from stellar evolution in the ini- 
tial phase of the H-ingestion phase for a model like those in 
(Herwig 2001), recalculated with f v = 30 and higher time res- 



olution. The proton abundance at any location is the result of 
mixing and simultaneous burning. The two times correspond to 
panel A and B in Fig. 4 in Miller Bertolami et al. (2006 ) and the 
account of events given in their Sect. 3.2.1 applies here as well. 



At time to the He-shell flash convection zone is about to 
make contact with the H-rich layers above. The H-profile at 
m r ~ 0.6042 M Q is the burning profile of the now extinct H- 
shell. During the late phase of the post-AGB evolution, ba- 
sically past the 'knee' in the HRD, the H-shell is inactive, and 
the He-shell convection can grow into the H-rich layers and mix 
those protons (and 3 He) down into the 12 C-rich He-shell flash 
convection zone. As H is mixed into deeper and hotter regions 
its lifetime against capture by 12 C decreases because the rate of 
the nuclear reaction 12 C(p,7) 13 N increases strongly with tem- 
perature. At some depth, in our simulation at m T = 0.6005 M Q , 
the mixing time scale equals the nuclear time scale (Damkohler 
number Da ~ 1, Sect. j 1 . 1 [ > and protons are now reacting rapidly 
with 12 C, thereby releasing for a brief period more energy than 
the He-shell that is intially driving the flash. 

In the stellar evolution simulation we treat time-dependent 
mixing mathematically as a diffusion process. It is implicitly 
assumed that on spheres the H-abundance is exactly homoge- 
neous, and that the radial mixing efficiency based on the radial 
mean convective velocity is also exactly homogeneous. This 
assumption in combination with the strong temperature sensi- 
tivity of the p-capture reaction causes the stellar evolution code 
to predict the shell of peak H-burning energy release to be ex- 
tremely thin. In the stellar evolution code an entropy step de- 
velops that separates the H-ingestion top convection from the 
He-shell flash convection underneath. A thin radiative zone 
formally prohibits mixing between the two convection zones. 
It shows up as a break in the diffusion coefficient line for time 
t\ in the top panel of Fig. [2] It now depends on the convec- 
tive boundary mixing assumptions whether or not material from 
the top convection zone can mix below and vice versa. These 
boundary mixing assumptions, i.e. the amount of overshooting 
appropriate for this situation, is not yet known. 

Fig. [2] shows that the split of the two convection zones ap- 
pears already very early when only a small amount of protons 
has been consumed. We mark the position in the H-profile 
and the corresponding H-abundance that has been reached at 
the time when the split occurs in the lower panel. The good 
agreement of our evolution simulation with the result by Miller 
Bertolam i et al.| ( |2006| Fig. 3 in their work) only means that 
these calculations properly converge and are precise, but not 
that they are accurate. 

At the time of the split the peak temperature in the now sepa- 
rated top H-burning driven convection zone is T < 1 .0 x 10 8 K. 
Although the 12 C(p,7) 13 N(/3 + ) 13 C reaction chain is providing 
plenty of the neutron source isotope 13 C, the 13 C(a,n) 16 reac- 
tion activation depends on the peak temperature reached in this 
top convection layer. For T = 10 8 K the lifetime of I3 C against 
capture by 4 He (and thus the time-scale of releasing neutrons) is 
454 yr, and thus neutron capture nucleosynthesis is negligible, 
considering that the born-again life time is only a few years. As 
a result, these stellar evolution models cannot provide the en- 



vironment to generate abundance patterns as observed by As 



plund et al. ( 1999 1 



We have performed a full nucleosynthesis analysis of the 
stellar evolution model sequence shown in Fig. [2] using the 



MPPNP code (Sect. A.l i. The technique for this nucleosynthe- 
sis analysis is explained in full detail in Sect.[5] Indeed no mod- 
ification of heavy element abundances is seen, in disagreement 
with the observations by Asplund et al. ([1999), and in agree- 
ment with the qualitative arguments that these authors made in 
their original paper. 
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The Herwig et al. (1999) models show a larger peak- 
temperature of T = 1.5 x 10 s Kp]for the H-ingestion driven top 
convection zone. As discussed in detail in Herwig (2001 ), those 
older models are not correctly reproducing the fast luminosity 
rise time observed in Sakurai's object, and it exists an inverse 
correlation between the rise time and the depth of the burning 
zone and split (i.e. convection speed, peak temperature), mod- 
els with the higher peak temperature have far too slow rise times 
and can thus be excluded. For these higher peak temperatures 
the life time of 13 C is 0.13yr. However, even this is not short 
enough to generate the abundance patterns observed in Saku- 
rai's object (see Sect.[5]for further discussion). 

We conclude from this analysis that a one-dimensional stel- 
lar evolution calculation cannot fully account for the mixing 
conditions in the convective-reactive H-ingestion flash that oc- 
cured in Sakurai's object. In this section we have already hinted 
at the possible reasons for the decrepancy. We will now have 
a closer look at what information and guidance we can derive 
from present three-dimensional hydrodynamic simulations of 
He-shell flash convection. 

4. THE HYDRODYNAMIC PICTURE 
4. 1 . New simulations 

In order to study the hydrodynamic process of entrainment 
and further mixing of H-rich material from the stable layers 
into the convection zone we have carried out new gas dynamics 
simulations of the entire three-dimensional He-shell flash con- 
vection domain in An geometry (Fig. [3}. We used the PPM code 
described in Sect. A. 2 We have not included burning of protons 
with 12 C because we restrict the goal of the numerical experi- 
ments purely to the investigation of mixing properties during 
the onset of the H-ingestion, which starts when the He-shell 
flash convecion has reached its largest Lagrangian extension. 

|Herwig et al.| f2006) simulated the He-shell flash convection 
shell as plane-parallel box-in-a-star. They selected an earlier 
phase of the He-shell flash when the convection had not yet 
reached its largest extent, and the H-rich layers had not been 
reached. Therefore only ~ 4.5 pressure scale heights needed 
to be included in those simulations which made them consider- 
ably less demanding than the new simulations. In addition the 
previous simulations were only in 2D. 

The new simulations were performed on a cubical domain 
with two uniform Cartesian grids of 576 3 and 384 3 respectively 
(Fig.[3]l[^] Each simulation realistically represents the abun- 
dance mixture in the He-shell flash convection zone and in the 
stable layer above as different materials with the correct molec- 
ular weight ratio. The setup includes an inert white-dwarf-like 
core and a radiative region below the bottom of the He-shell 
flash convection zone at 9, 500km where the gravitational accel- 
eration is 4.9545 ■ 10 7 cm/s 2 , the density is 1.174- 10 4 g/cm 3 and 
the pressure is 1.696 ■ 10 g/cms 2 . At the bottom of the convec- 
tion zone a luminosity of 4.2 • 10 7 L Q is artificially added in a 
shell of 1 , 000km. This heating corresponds to the He-burning 
that drives the flash, and compares as follows to the He-shell 
flash luminosity in the stellar evolution models shown in Fig. [2] 
In the model at time fo the He-burning luminosity is at its peak 

1 1 We have now recalculated those old models with higher resolution and find the peak H-burning location at slightly lower temperature of T = 1 .3 X 10 s K. 

12 The 576 3 calculation took 4 days on 24 workstations at th e University of M innesota's Laboratory for Computational Science & Engineering (LCSE). A movie 
made from the output of this run may be downloaded from the |LCSE Web site] 

Ll We have continued this run for another 14 convective turnovers. However, as will become clear from the following discussion the omission of proton burning limits 
the scientific use of that later part of the run to our application. Note that the time step of the 3D simulations is limited to At = 5.9 ■ 10~ 2 s which implies that 300,000 
cycles had to be computed to reach the state shown. 



of Lee.o = 4.75 • 10 7 L© whereas it drops somewhat once the H- 
burning flash ignites at t\ when L He ,i = 4.27 • 10 7 L Q . Thus, the 
3D hydrodynamic simulations are driven at the nominal heating 
rate. 

The top of the convection zone is at a radius of 30,000km 
and surrounded by a radiative shell of thickness 4.500km. The 
three layers are each polytropes. The adiabatic polytrope that 
represents the convection zone spans ~ 9H p . The setup con- 
tains two materials. The lighter material represents the H/He 
mixture in the stable layer above the convection zone. The 
heavier fluid represents the I2 C-rich mixture that occupies the 
convection zone. We have assumed here that the material in the 
stable layer below the convection zone has the same molecular 
weight. The ratio of the molecular weights of the two compo- 
nents is ^c,0,He/A*H,He = 2.26. 

The higher resolution run (Fig. [3] right panel) is shown at 
time 21,653s. For convective transport the typical radial ve- 
locities are of interest. In the shown snapshot the largest ra- 
dially rms-velocities are found about 4, 500km above the bot- 
tom of the convection zone around < v ra d. aV e >= \/2 < £km > ~ 
12.5km/s. The velocity of individual convective gusts can be 
significantly higher. Towards the upper boundary of the con- 
vection zone the radial velocities decrease to a few km/s. This 
is compensated by large tangential velocities > 12km/s which 
stay this high all the way to the convection boundary (Fig.Q. 
The resulting strong radial gradient of the tangential velocities 
at the top convection boundary is, via Kelvin-Helmholtz in- 
stabilities, likely the main mechanism of the entrainment and 
convective boundary mixing that we observe in these simula- 
tions. The information on typical convective velocities together 
with the radial scale of the convection zone implies a convec- 
tive turn-over time scale of the order ~ 3000s (cf. Appendix 15), 
Therefore, Fig. [3] shows the entrainment after ~ 7 convective 
turnovers When estimating the time scale for H-rich mate- 
rial to enter the convection zone it must be considered that the 
entrained material is dominantly transported in downflow lanes 
that are gravitationally compressed as the material descends. 
This mechanism is reflected in the radial velocities of the H- 
rich material that has entered the convection zone, which in the 
snapshot shown exceed 20km/ s. We note that for this compo- 
nent even the radially averaged velocity corresponds to a Mach 
number Ma ~ 0.02 which is much higher than the MLT con- 
vective velocity based estimate of Ma ~ 0.001. 

After some initial transient period the convection assumes 
a flow pattern that is dominated by large upwelling convec- 
tive cells that occupy typically a full octant as they emerge at 
the top convection boundary. These large convective structures 
can be observed because we simulate the full 4tt sphere. En- 
trainment of the H-rich material from the stable layer into the 
convection zone is mostly associated with downdraft lanes that 
form when large cells collide on the surface of the convection 
zone (Fig. [3]). Note that the radially averaged profile of the in- 
gested H-rich material from the 3D hydro-simulation is qual- 
itatively very similar compared to the diffusion picture of the 
one-dimensional stellar evolution (Pig.[5j, at least close to the 
upper boundary. Futher inward the lines divert from each other 
systematically as no H is burned in the 3D simulations (this 
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physics is not yet included). 

However, the important result of the 3D simulations is that 
entrainment is rather inhomogeneous and asymmetric, as well 
as intermittent in locally confined wedges of the star. From 
the snapshot image of the entrainment it is clear that signifi- 
cant anisotropy of the H-abundance is advected into the deeper 
layers where the burning will eventually take place. 

4.2. Implications for the nucleosynthesis in a convective 
reactive environment like Sakurai's object 

We will give a full account of these simulations elsewhere. 
Here we want to describe a few properties that are relevant for 
guiding our mixing strategy for the nucleosynthesis simulation 
of the flash in Sakurai's object. The details of the convective- 
reactive burning of hydrogen in the He-shell flash convection 
zone depend on two aspects of the problem that hydrodynamic 
simulations can address. The first is the process of entrainment. 
How much is the fuel is premixed immediately after the en- 
trainment in the near-boundary layers. Subsequently these H- 
enriched fluid elements will be carried along with the convec- 
tive flow to deeper and hotter layers where protons will even- 
tually react with 12 C. This leads to the second aspect of the 
problem, the hydrodynamic feedback of the nuclear energy re- 
leased. In the one-dimensional simulations this feedback is in 
the form of a sharp entropy barrier, or a thin shell of positive en- 
tropy gradient locally confined to a sphere. In reality the thick- 
ness of this layer will depend on the velocity distribution and 
the abundance distribution of fluid elements entering the layers 
hot enough for rapid burning. 

We can illustrate the possible outcomes by considering two 
extreme cases. Assuming first that any entrained material is im- 
mediately mixed and that vertical velocities of fluid elements 
are only deviating negligibly from some average value (obvi- 
ously, this case is very close to the MLT picture of convection) 
then all fluid parcels or blobs would release nuclear energy at 
almost the same radial position inside the convection zone, and 
thus a very thin burn layer would form, concentrating the en- 
tropy jump into a narrow region with large positive entropy 
gradient, and soon inhibiting any further radial mixing. The 
other extreme would be a wide range of mixing ratios in blobs 
of H-enriched material entering the deeper layers with a large 
range of velocities. Both of these inhomogeneities lead to a 
broadening of the burning layer. To first approximation a blob 
(note that this may be a shredded blob in order to conceptually 
overcome mixing-length concepts) burns at Da ~ 1 (Sect. fLjj , 
For smaller Da (above the burning layer) the nuclear reaction 
time scale is longer than the mixing time scale and the blob 
will rather move further down than burn. For Da > 1 we are 
below the burning layer because now the blob burns faster than 
it can move further down. Since the burn time scale decreases 
with depth a range of mixing velocities translates into a spatial 
range in which Da ~ 1. Differently than in the first case, the 
velocity distribution of blobs leads to a broadening of the burn 
layer. Distributing the energy released from proton capture over 
a thicker layer will make the emerging entropy gradient shal- 
lower. Mixing accross the burn layer will be more efficient. A 
distribution of levels of H-enrichments in blobs being advected 
through the burn layer would mean that the H-abundance is het- 
erogeneous (patchy) on spheres. Thus, the energy generation 
and the dynamic feedback may very well be patchy and inho- 
mogeneous on spheres, as well as time variable. At least ini- 
tially, the inhibiting effect of the burn layer on mixing may as 



well be time variable and inhomogeneous on spheres. 

In other words, an inhomogeneous distribution of fuel abun- 
dance in blobs together with a distribution of vertical blob ve- 
locities would have the tendency to delay the inhibiting effect 
of nuclear burning on mixing from the top to the bottom of the 
convection zone. We leave a detailed quantitative analysis of 
these processes to a forthcoming investigation. Here we focus 
on the conceptual guidance we can gain from the hydrodynamic 
simulations. These do indeed show a significant inhomogeneity 
of the entrained material all the way down to the bottom of the 
convection zone (Fig. [3}, as well as a significant distribution of 
vertical velocities, including convective gusts up to Mach num- 
bers around Ma ~ 0.03. 

We conclude from this analysis that the hydrodynamic na- 
ture of the convective-reactive phase of H-ingestion into the 
He-shell flash convection zone likely translates into a contin- 
ued mixing through the burn layer. We therefore hypothesize 
that mixing is not inhibited at the early stage, as indicated by 
stellar evolution models, but that instead mixing accross the H- 
buming layer is possible for a prolonged period. It may stop 
only at a later time after more H-ingestion has taken place. In 
the next section we will test this hypothesis through nucleosyn- 
thesis simulations that can be compared with the observations 
Asplund et aTT] (| 1 999|> . 



5. NUCLEOSYNTHESIS SIMULATIONS 

In this section we will describe mixing and nucleosynthe- 
sis simulations based on the thermodynamic stellar evolution 
structure of a post-AGB He-shell flash. We describe intially 
two cases, one that resembles the mixing predicted by stellar 
evolution (Sect. 5.2 1, and one with a mixing prescription that re- 
flects the findings discussed in the previous section (Sect. 5.3 i. 
While the first fails to reproduce key observational features of 
Sakurai's object, the second one succeeds. We show that high 
neutron densitities in the range 10 12 < Af n /cm" 3 < 10 16 are re- 
quired to reproduce the observed abundance features, as already 

T999J. Such a neutron density 



pointed out by Asplund et al. 



regime is higher than the classic s process and significantly 
lower compared to the classic r process. 

5. 1. General setup of nucleosynthesis simulations 

We are using the MPPNP post-processing code (Sect. |A.l"| 
to calculate the nucleosynthesis of a He-shell flash peak one- 
dimensional stellar structure model. We use two structures, one 
of them shown in Fig.[2]for t = tQ. The MPPNP code reads the 
mixing-length theory diffusion coefficient as well as the tem- 
perature and density structure from the stellar evolution struc- 
ture model. We post-process this structure with sub-time steps 
of Afpost-processing = 63s. Thus, the mixing time scale is well re- 
solved, and the numerical splitting of the mixing and the nucle- 
osynthesis operators is justified. The He-shell flash convection 
zone is spatially resolved with 70 to 90 zones. The grid is stat- 
ically refined and provides extra resolution near the ingestion 
layer at the top of the convection zone, as well as around any 
split region, should it be included. 

The MLT based diffusion coefficient that is read in along 
with the stellar structure from the stellar evolution output does 
not show a split because the stellar evolution model is from a 
time just before the ingestion of H-rich material begins. How- 
ever, we are providing for an optional split that can be inserted 
at an arbitrary location and time, by modifying the diffusion 
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coefficient in Eulerian coordinates in the following way 

Dmlt 

(1. +a 2 exp(-a i (m r - m r , sp iit) 2 ) 



D 



with split " 



(2) 



where the split is located at ra, -, sp ii t . When a split is imposed it is 
chosen to be deep enough so that only very little material can be 
mixed through, and the split is also very narrow. With a\ = 10 4 
and ai = 10 7 the diffusion coefficient in the convection zone 
of £>mlt ~ 5 • 10 13 cm 2 /s is reduced to D sp n Um \ n ~ 5 • 10 6 cm 2 /s 
over a width of < 10~ 4 M Q . We emphasise that a\ and are 
free parameters of our simple delayed split model and their par- 
ticular value is not important at this point. Only further hydro- 
dynamic simulations can possibly determine the mixing prop- 
erties in this environment. The purpose of the delayed split in 
terms of the radially averaged nucleosynthesis calculations is 
further discussed below. 

We are solving only for the nucleosynthesis and mixing equa- 
tions while the T, p stratification is assumed to remain un- 
changed. Protons and 3 He are inserted into the top of the con- 
vection zone at a rate that is derived from the Lagrangian veloc- 
ity of the top of the convective boundary, as it moves into the 
H-rich layers above the convection zone in the stellar evolution 
model. This velocity is M to p,conv ~1.7x 10~ 2 M Q /yr. We are 
ingesting at a rate of 5.3 x lO" lo M /s[^] We also add 3 He ac- 
cording to the solar H/ 3 He ratio in order to obtain a prediction 
for Li. 

Another constraint is that the total amount of H available for 
ingestion is limited to the small remaining envelope mass that 
remains on the pre-formed WD when the star leaves the AGB. 
For a core mass of O.6M this envelope mass is ~ 1O~ 4 M 
with H and He fractions as expected at the end of the AGB 
(mostly the initial ratio possibly modified by third dredge-up). 
In all of the cases discussed here we always find a nucleosyn- 
thetic reason to stop a simulation before we run out of fuel. 

5.2. Stellar evolution mixing case 

In the stellar evolution models the convection zone split 
due toH-burning activation starts as soon as H is ingested 
(Sect. 3.2 1, and no H or 13 C can by mixed below the split co- 
ordinate. In Fig. [6] we show the abundance distribution pre- 
diction at the top of the convection zone for this model in com- 
parison with the observations by |Asplund et al.| ( |1999[ l. We 
have used the (p,T,D) stratification (strat-A) from the ( |Her-| 
wig et al. 1999| > sequence, selecting a model just before the 



flash convection zone from 22 Ne(a,n) 25 Mg reaction. Neverthe- 
less, the predicted abundances do not matching the observa- 
tions. 

Li was produced initially during the ingestion (see below and 
Herwig & Langer 2001 for more details) is destroyed on the 



time scale of ~ 1 yr. Stellar models predict that material around 
and beyond the split expands and cools which reduces the a- 
capture efficiency depleting Li. But this also reduces the pro- 
duction of heavy elements. 

Sc is well reproduced within the uncertainties, in neutron 
densities higher than in the classic s-process. 40 Ca is the main 
seed along the neutron capture path, and Sc is mainly produced 
as 45 Ca which will decay to 45 Sc in ~ 166 days. The production 
of Sc is subject to nuclear reaction uncertainties, for instance 
the (n,7) rates of Ca isotopes, 41 Ca(n,p) 41 K and in particular 
4I Ca(n,a) 38 Ar. 

The bottom line is that Li observations cannot be repro- 
duced together with a significant s-process nucleosynthesis in 
this simulation. But most importantly, the predicted [hs/ls] ra- 
tio much higher than observed. Therefore, the nucleosynthesis 
simulation based on the one-dimensional stellar evolution pre- 
diction for mixing cannot account for the observed abundance 
patterns in Sakurai's object, which confirms our findings from 
Sect.l3~2l 



5.2 



4.2 1. We use same background model (strat-A) as in 



H-ingestion starts as a template for this run. The mixing split 
as described in the previous section is activated immediately as 
H starts to mix into the convection zone. Peak H-burning is 
located at a higher temperature in the ( |Herwig et al.| 1999] > se- 
quence compared to more recent models, and therefore this case 
yields an upper limit of the nucleosynthesis efficiency predicted 
from one-dimensional models. 

Calculations are run for about one year, after which also Ba 
starts to be produced, in disagreement with observations. The 
neutron density reaches a value of the order of 10 11 cm" 3 at 
the split coordinate due to the high 13 C concentration accumu- 
lated via proton capture on 12 C. This value is comparable with 
the neutron density obtained at the bottom of a regular He-shell 

14 Specifically, we add every ~ 6min (every 6 lh cycle, corresponding roughly to 10 times per convective turn-over time) AX = 5 X 10~ 4 to the mass fraction of H in 
the uppermost 4 X 1O~ 4 M0 of the con vectio n zone. The baryon numbers are conserved by subtracting the required mass fraction from 12 C. The abundances up to 
23 Na are initialized as described in Sect. 3. 1 

predict 



5.3. Delayed split model motivated by the hydrodynamic 
simulations 

We now assume that the split is not created instantaneously 
by H-burning, but mixing continues — at least initially — un- 
restricted despite the energy generation from H-burning (see 
Sect. 

Sect, 

13 N is still formed in the upper layers where the reaction and 
the mixing time coincide (Fig.[7jl. 13 N decays to 13 C on a time 
scale of ~ lOmin. During this time 13 N will be swept along with 
the flow, possibly covering a distance of the order 10,000km. 
Eventually 13 C is mixed to the bottom of the He-shell flash con- 
vection zone (T ~ 2.5-3.0 • 10 8 K) and establishes an abun- 
dance of ~ 1% by mass throughout the He intershell. Neutrons 
are released via 13 C(a,n) ie O on the time scale of 1 . . . 10s and 
neutron densities reach a value of ~ 10 15 cm at the bottom of the 
convection zone. The profile for Sr is shown as an example for 
how the abundance, even of heavy elements, varies inside the 
convection zone as mixing and production proceed at similar 
time scales. 

The intense neutron burst leads to the formation of the first s- 
process peak elements Rb, Y, Sr, Zr, with Fe as the main seed. 
The unimpeded mixing between the formation region of 13 N 
and the deeper layers where the neutrons are released must fin- 
ish before the Ba-La elements are significantly produced, which 
is not observed. This defines the moment when mixing finally 
has to be limited, and we then turn on the delayed split. In 
Fig. [8] we show the abundances expected at the top of the He 
intershell for different split time between 800min to 1200min. 
Burning of 3 He produces 7 Be via the reaction 3 H e(q,7) 7 Be, 



which will decay later to 7 Li. As pointed out by Herwig & 
Langer (2001 ), Li destruction is avoided under these conditions 
not because Li is mixed into cooler regions (Cameron-Fowler 

-4, 



15 



As discussed in Sect 



1.2 



this older model did not reproduce the observed light curve, but more recent models Herwig 



the split at lower temperature and as a result even less n-induced nucleosynthesis. 



(2001 



Miller Bertolami et al. 



2006 
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mechansim). Rather, in this hot H-deficient 3 He-burning all 
the protons are consumed before 7 Be decays to 7 Li. Then 7 Li 
is more stable as it is only destroyed through a-captures. In 
all cases Li is overproduced if we can assume that a sufficient 
supply of 3 He is still available in the envelope when the VLTP 



begins (cf. Sect. 1.2 1 



Mg is more abundant in the simulations by one order of mag- 
nitude compared to observations. In all runs Mg is only weakly 
modified by nucleosynthesis. For this reason, the low observed 
Mg abundance may be another indicator of a sub-solar initial 
metallicity of the star, unless there is some observational prob- 
lem. 

Despite the differences between these tests and the measure- 
ments, the overall abundance trends are similar. In particular all 
three test cases in Fig [8] have a low [hs/ls] that ranges between 
-0.9 and -1.5, decreasing with increasing the split delay. In the 
case with the latest split (at 1200min), [hs/ls] is still ~ 0.5dex 
higher than observed in Sakurai's object. Light and intermedi- 
ate elements are not much affected by the split time. 

In addition to the split delay time the quantitative model pre- 
dictions depend on the base stratification and convective mix- 
ing coefficient taken from the stellar evolution model. This de- 
termines, for instance, how quickly the protons and resulting 
13 C are mixed, and in turn the neutron density. To test the de- 
pendence of the results on this point we present another set of 
simulations based on the structure (strat-B) at the last thermal 
pulse in the 2M Q star model sequence by Herwig et al. (model 



ET14 2006). We have applied a delayed split as for the strat-A 
model. With this base structure, the measured [hs/ls] is repro- 
duced within the uncertainties (Fig. [9). However, now Zr is 
higher by 1 dex compared to the Asplund et al. measurements. 
A general overview of the abundance profiles in the He inter- 
shell for the most indicative light isotopes and of the elements 
included in Fig. [9]is given in Fig. 10 where the split position 
and the variation in the abundances are shown. 

To] (left upper panel) confirms that the 12 C/ 13 C = 6.7 



Fig. 



ratio agrees within uncertainties with the observed ratio of 
~ 3 ... 5. 7 Be is shown in the lower left panel to be highly abun- 
dant, which will feed Li. 

This profile view of one of our simulations reveals that 
the neutron capture nucleosynthesis continues below the split, 
thereby further modifying chemical abundances. Possibly this 
further processed material below the split has affected Sakurai's 
observed surface abundances, through additional, later mixing. 
H-burning at the split must lose efficiency at some point when 
running out of fuel. This may allow material exchange between 
the two regions (see also discussion in Asplund et al. 1999 ). As- 
plund et al. observed Sakurai's object four different times in 6 
months, and these observations show some drastic changes for 
some elements. It is not the aim of this paper to directly address 
these abundance trends over the sixth-month period, since this 
level of detail cannot be captured by our modeling approach, 
but has to await updated multi-dimensional simulations. 

However, we may assume, as a working hypothesis, that the 
He intershell is made of two components, one heavily processed 
below the split (region 1), and one above the split (region 2) 
that was affected only by the first ingestion phase. Because of 
the decreasing of efficiency of the H-burning at the split, some 
material from region 1 is allowed to reach region 2 again and 
contribute to the observed abundance distribution. Such a two- 
component model is shown in Fig.[TT] Starting from the simu- 
lation based on stratification strat-A, with a delayed split after 



1200min (see Fig. [8} 10% of the material is coming from re- 
gion 1, and 90% from region 2. No significant differences are 
obtained compared to Fig. [8] However, this depends on how 
much material is mixed from region 1 to region 2. In this spe- 
cific case, such mixing implies a decrease on [hs/ls], but also 
an increase on Ba production, not supported from the observa- 
tions. For this reason, at present we cannot confirm or rule out 
such a double component scenario. 

5.4. Nuclear reaction rate uncertainties 

In this nucleosynthesis scenario both H- and He-burning re- 
actions, as well as the n-capture reactions including those of 
short-lived isotopes, are important. Especially, several elemen- 
tal abundances, for example Ti and Sc, are strongly dependent 
on s-process branchings which requires extra accuracy from the 
nuclear physics data. As we want to use this case to probe fu- 
ture hydroynamic simulations we need to asses the influence of 
nuclear rate uncertainties. 

In Fig. [12] we show for the model strat-A with split af- 
ter lOOOmin the impact of changing the 13 C(a,n) ie O and the 
14 N(n,p) 14 C reactions by a factor of 2. The 25 Mg(n,7) 26 Mg 
reaction has been varied by a factor of 1.2. 13 C(a,n) 16 is 
the main neutron source and the two neutron capture reactions 
are important neutron poisons. Among these tests, the [hs/ls] 
changes between -0.9 and -1.6. In particular the first peak 
elements are strongly affected. The Rb abundance changes 
by ldex. Intermediate and light element predictions are only 
weakly affected by nuclear reaction rate uncertainties. Small 
errors associated with the CNO cycle rates (e.g., 12 C(p,7) 13 N 
and 14 N(p,7) 15 0) have a marginal impact in our results com- 
pared to the other rates that we have considered. 

In Fig. 12 we only included the impact of varying the neu- 



tron capture reaction rates of light neutron poisons. In the short 
time scale of the neutron burst, the neutron capture process 
is also expected to show a strong propagation effect in the fi- 
nal abundance distribution beyond iron, due to uncertainties of 
neutron capture rates along the nucleosynthesis path. In partic- 
ular, such propagation may be relevant in our case, since Rb, 
Sr, Y and Zr production is affected by the error of several low 
cross sections of isotopes in the mass region between Fe and 
Sr, acting like bottle-neck s in the neutron capt ure flow (e.g., 
62 Ni, 68 Zn, 74 Ge and 78 Se |Pignatari et al.|2010] and reference 
therein). Another point to consider is that in the high neutron 
density regime reached in our calculations several unstable iso- 
topes are produced efficiently, and many stable isotopes receive 
a significant contribution from unstable species from radiogenic 
decay and/or from decay during the neutron freezout, when the 
split is established. For instance, in all the cases presented in 
Fig. [12] most of Y (that is formed by one stable isotope only, 
89 Y) is produced as 89 Sr. The neutron capture rates of unstable 
species are mostly theoretical, and also their large uncertainty 
(typically a factor of 2-3) may affect the final isotopic distribu- 
tion. 

None of our simulations seem to be reproducing Sc partic- 
ularly well. Sc and the elemental ratio Sc/Ca are particularly 
sensitive to the neutron density. Indeed, 45 Sc is produced as 
unstable 45 Ca via neutron captures on stable Ca species, where 
40 Ca is the main seed for Sc production. 41 Ca is unstable, and 
has stronger (n,p) and (n,a) than (n,7) channels. For this rea- 
son, the uncertainty in the relative efficiency of the (n,p), (n,a) 
and (n,7) channels may affect the total Sc production. Among 
nuclear uncertainties, another possible explanation for Sc over- 
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production is that the initial metallicity of the Sakurai's ob- 
ject is even lower than what we have used for our simulations 
([Fe/H]=-0.18). Indeed, a lower initial 40 Ca will results in a 
lower final Sc abundance. 

6. CONCLUSIONS 
6.1. Summary 

We have presented in this paper a multi-physics view of 
the combustion in a very-late thermal pulse in a pre-WD. H 
is mixed convectively into the He-shell flash convection zone. 
We have discussed the one-dimensional stellar evolution pic- 
ture, that predicts that early on the energy generation from 
the 12 C(p,7) 13 N reaction creates a sharp entropy discontinuity 
which prohibits mixing. A detailed nucleosynthesis analysis, 
based on a complete multi-zone treatment of nucleosynthesis 
with mixing, shows that this one-dimensional structure evolu- 
tion leads to abundance predictions that are incompatible with 
the observed abundances in Sakurai's object. Seeking guidance 
from full three-dimensional hydrodynamic simulations of He- 
shell flash convection in 4ir geometry with entrainment, we ob- 
tain reasons to suspect that the burning front is more distributed 
than predicted in one-dimensional stellar evolution. Fuel will 
be transported down in down-draft lanes leading to an inhomo- 
geneous distribution of fuel in the burning zone. In addition, 
vertical down drafts enriched with fuel will populate a velocity 
distribution. From this information we speculate that mixing 
of protons and of the neutron source material 13 N which later 
becomes 13 C, across the convective H-burning zone will pro- 
ceed for much longer than indicated by one-dimensional stellar 
evolution. 

We point out that the main nucleosynthetic signature of 
convective-reactive burning in this study is the significant over- 
production of the first peak elements Sr, Y and Zr, coupled with 
a non-efficient production of the second peak elements Ba and 
La. According to our simulations, neutron densities 10 12 cm" 3 
< N n < 10 16 cm" 3 are required to explain such abundance dis- 
tribution. More specifically, in the Sakurai's object time scale 
of ~ 2 years between the luminosity peak due to H burning 
and the Asplund's observations, a neutron density peak of ~ 
10 15 cm" 3 with a delay of ~ 1 day before the complete split ac- 
tivation would qualitatively reproduce the observed [hs/ls] ra- 
tio, the Li abundance and the low 12 C/ 13 C ratio. The problems 
that we encounter in reproducing single elements may be due 
to the approximations in our model (e.g., for the nucleosynthe- 
sis simulations we use parameters from one-dimensional stellar 
models), to observation problems (e.g., the observed Y/Zr ratio 
cannot be reproduced by neutron capture nucleosynthesis) or to 
nuclear physics uncertainties (e.g., Sc). 

Nuclear reaction rate uncertainties are shown to have a par- 
ticularly important effect on some key observables in this non- 



equilibrium nuclear burning environment. 

6.2. Implications for stellar evolution and nucleosynthesis of 
the first generations of stars 

One of our main motivations to study convective-reactive 
phases in stellar evolution is their prevalence in models of the 
first generation of stars. As reviewed in Sect. jl . 1| convec- 
tive mixing of protons with the I2 C from He-burning at He- 
buming temperatures is frequently encountered in stellar evo- 
lution calculations at very low and zero metal content at all 
masses. This investigation shows that the predictive power of 
one-dimensional stellar evolution simulations is severely lim- 
ited for observables that depend on these convective-reactive 
phases. 

Neutron burst nucleosynthesis of the type described in 
this paper are nevertheless expected to happen also in the 
convective-reactive H- 12 C combustion events in first generation 
of stars. Indeed, the neutron source 13 C is of primary origin, 
i.e. its abundance is not affected by the metal content in the 
initial stellar composition. Massive stars at different metallic- 
ities may experience H- I2 C combustion, ingesting protons in 
the He shell (see discussion in Woosley & Weaver 1995|. If 
enough hydrogen fuel is ingested then Sr, Y and Zr may be ef- 
ficiently produced by the primary 13 C(a,n) 16 neutron source, 
just as in our simulations presented here. This may be an alter- 
native or complementary explanation for a missing component 
in the first neutron-peak region of the abundance distribution in 
both the solar abundance distribution as well as the metal poor 



stars, (light-element primary process, or LEPP 


Travaglio et al. 


|2004||Pignatari et al.|2008||Farouqi et al.[2009 


i. In the future 



we intend to study the speculation that the convective-reactive 
proton- 12 C combustion in the convective He shell in massive 
stars could provide another possible solution for the LEPP. 
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FIG. 1 . — Observed and predicted s-process abundance distribution index ratio [hs/ls] for stars with a large range of metallicities. Observa tions I Tech 1971 ; Smith 
1984, Smith & Lambert 1984 1985. 1986. Smith & Suntzeff 1987, Smith & Lambert 1990, Smith et al. 1993 1996 1997, Abia & Wallerstein 1998, Van Winckel 
& Reyniers 2000 Zacs et al.|1 995 1998 , Zacset al.|20 00 Reddy et al. 1999 Kipper et al. 1996 Kipper & Jorgensen 199 4||Tomkin & Lambert|1 983 1986 Kovacs 
1985|| Vanture|1992||2000||Pereira et al.|1998||Aoki et al.|2U 00 McWilli arnet al.|1995)|McWilliam|1998||Norris et al.|1997||Beveridge & Sneden|1994) and model 
predictions of AGB stars are from Busso et al. (2001 1. In the Figure, the [hs/ls] ratio observations of the Sakurai's object have a certain range, depending on which 
of the 4 observations from Asplund et al. are considered, and how the indices are calculated. In general, the ratio is about 2dex smaller compared to AGB predictions 
and observations. Our nucleosynthesis results are also included for comparison (see Sect.[5]for details). 



Table 1 

Observed neutron capture signature, |Asplund et al.| i |1999| >. 



[Fe/H] = (-0.63) 


April 1996 October 1996 


[Y/Fe] 


+0.96 (+1.59) +1.96 (+2.59) 


[Ba/Fe] 


-0.63 (0.0) -0.23 (+0.40) 


[Ba/Y] 


-1.59 (-1.59) -2.19 (-2.19) 
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FIG. 3. — Hydrodynamic picture of H-entrainment into He-shell flash convection near the luminosity peak of the flash. The setup is based on a stellar evolution 
model corresponding to the situation shortly after time shown in Fig.[2] when the top of the convection zone is just making contact with the H-rich stable layer. 
Colors indicate abundance of proton-rich material that is originally only in the stable layer above the convection zone that is entrained into the convection zone. 
Volume fractions of about ~ 1 % are shown as blue, while concentrations that are close to one are transparent. The lowest concentration yellow blobs that are mixed 
deep into the convection zone correspond to ~ 0.01%. Abundance levels below approximately 5 X 10~ 5 have been made transparent as well. The left panel shows a 
snapshot from a 384 3 grid while the right panel image is from a ru n on a 576 3 grid. Slightly different times are sho wn and similar but not identical color maps have 
been used. The PPM simulation is described in more detail in Sect. |4.1| and the simulation code is described in Sect. |A.2| 
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FIG. 5 . — Comparison of entrainment of material from the stable layer above the convection zone into the 12 C-rich layer as it is represented in the one-dimensional 
stellar evolution model with mixing treated as diffusion in the mixing-length picture and in the 3D simulations discussed in this paper. The 3D profile (green line) 
shows the same data, radially averaged, as in Fig.[5] right panel. The ID line (blue) is the line labeled to in Fig.[2] The mass coordinates have been set to zero in both 
cases near the top of the convection zone. 
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FIG. 7. — The abundance profiles snapshot (RUN103), just before the mixing split is imposed, demonstrates the simultaneous action of nucleosynthesis and mixing 
on similar time scales. 
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Abundance distribution at step 2000 for different cases with the split starting after 800min (RUNI05), lOOOmin (RUN103) and 1200min (RUN106). 
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FIG. 9. — Abundance distribution at the end of the simulations (RUN48/strat-B) after 3000min, when all H- and 3 He-ingestion has been ingested. 
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further mixing between the He-shell flash driven convection zone (left) and the H-ingestion flash driven convection zone (right). Arrows in the right panels indicate 
the observed abundances, connecting the solar values with the observed ones. 
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FIG. 1 1 . — Abundance distribution at step 2000 for a split delay of 1200min, considering mixing of 10% of the deep component with 90% from the component 
above the split. 
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Abundance distribution for different nuclear test at step 2000, from RUN103 (split delay = lOOOmin) as standard, and RUN107,RUN108,RUN109. 
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APPENDIX 



CODE DESCRIPTION 
Nucleosynthesis 

The PPN physics package allows a flexible combination of nuclear reaction rates and entire compilations of rates. For this study 
we choose for the main charged particle reactions the compilation by Angulo et al. (1999) (NACRE compilation). This choice 
allows us to be consistent with the original network used to calculate the stellar stru ctures for basic energetic nuclear reactions, i.e., 

14* 



N(p,7) O, 3-a and C(a,j) O. Notice that the use of more recent rates (e.g., Imbriani et al. 2005 



Fynbo et al. |2005 Kunz 



et al. 2002 respectively) would not change our results, where uncertainties related to physics processes and mixing still has a critical 
impact. Other charged particle reactions, among the others 13 C(a,n) 16 0, which is the main neutron source during the H ingestion, 
have more recent measurements (e.g., |Heil et al.|2008] l. However, in this case NACRE rates are consistent with the new rates within 
their uncertainty. For instance, we consider a factor of two of uncertainty for the I3 C(a,n) 16 rate in the temperature regi me that is 
relevant for the 13 C burning (see Section[5]for more details). For neutron capture reactions of stable isotopes we refer to Dillmann 
et al. (2006) (KADoNIS compilation). Stellar /3-decay rates and electron captures are from |Odaetal.| ( [T994| and |Fuller etaL| ^ 1985) 



for many light unstable isotopes, and Goriely ( 1999 ) for many heavy unstable isotopes. Rates not included in the previous references 
are given by the Basel REACLIB compilation. We are solving the complete network in each radial grid point, including all relevant 
charged particle, n-capture reactions as well as the /3-decays. A recursive, dynamic network generation has been integrated into the 
solver, i.e. the size of the network automatically adapts to the conditions given. If, for example, a neutron source is activated the 
network will be automatically enlarged to include all heavy and unstable isotopes as needed according to the network fluxes. This 
dynamic network feature ensures that the network calculation never misses any important isotope or reaction. 

In these simulations we are using the multi-zone driver of the PPN code (MPPNP) that allows for the calculation of the complete 
nucleosynthesis in all of the zones of one-dimensional profiles, e.g. from stellar evolution, of density and temperature. The MPPNP 
driver employs MPI parallelism to enable efficient calculations on up to 30-50 processors depending on problem sizes. The simula- 
tions carried out here involve relatively small grids between 70 and 90 zones. A fully implicit nucleosynthesis step is followed by a 
mixing step according to the diffusion coefficient taken, for example, from the stellar evolution model. This procedure is repeated 
for subsequent time steps in order to compute the evolution of the abundance profiles of all species involved. Mixing and network 
calculations are performed in the operator split mode, which is a good approximation for the post-processing because we choose the 
post-processing time step to be small enough to resolve the mixing time scale. 

Hydrodynamics 

The PPM gas dynamics scheme (Woodward & Colella 1984 Colella & Woodward 1984 Woodward 1986 2006 ) has been in use in 
computational astrophysics for many years. It is incorporated in the community codes VH1 (jBlondin & Lufkin 1993 ), ENZO ([Bryan 
|et al.|1993] l, and FLASH ( |Calder et al.|20 02 ). The version that we use in this work is described in full in |Woodw ard (2006). Here we 
have augmented PPM with the PPB moment-conserving advection scheme to treat the entrainment of fluid from above the convection 
zone during the helium shell flash in an AGB star (see [Woodward et al.|2008| l. PPB is built upon van Leer's Scheme VI ( |van Leer] 
|1977] >, a 1-D scheme that conserves the first 3 moments of the advected distribution in each grid cell. To this scheme we have added 
a set of very carefully constructed constraints (Woodward 2005) keeping the advected fractional volume of a multifluid constituent 
of the gas within the range from to 1. These constraints are a considerable improvement over those outlined in Woodward 1986 for 
a 2-D PPB scheme. We have also streamlined the implementation of PPB in 3-D by eliminating various high-order terms in order to 
obtain a highly efficient, directionally split scheme (Woodward 2005 ) that conserves 10 moments of the distribution of the advected 
fractional volume variable in each cell. PPB is combined with PPM to describe multifluid hydrodynamics by adding the constraint 
of pressure and temperature equilibrium within each grid cell. At present our code is explicit. Mach numbers in the convective gusts 
of helium shell flash convection are about 1/30 or less. Consequently, we must take many time steps to follow the flow through an 
entire circuit of a large convective eddy. We note that such eddies are global in scale, and we follow them by including the entire 
convection shell in our computational domain. The conclusion that large scales are involved here is similar to the earlier findings of 
|Porter et aL] ( |2000| l and [Porter & Woodw ard (2006) for the outer convective envelope of an AGB star. The restricted time step values, 
from explicit hydrodynamics, and the large domain, arising from the natural scale of the convection, place significant demands on 
the computation. We address these demands in two ways. First, we exploit a new implementation of our codes aimed specifically at 
the multicore CPUs found in modern computers (see Woodward et al.|20 08 2009), which has delivered to our codes roughly a 40x 
speed-up over performance on single-core platforms from about 4 years ago (the code performance now stands at 24 Gflop/s/4-core- 
CPU, scalable to thousands of CPUs, and we obtain sustained performance over 1 Tflop/s on our small local cluster daily). Second, 
we exploit the fact that explicit computation is roughly as efficient as implicit computation when Mach numbers are around 1/30. 

The code scales to hundreds of thousands of processor cores, for which runs with the proper heating rates, the full convection 
zone, and well resolved entrainment at the convection zone boundary are easily carried out in a single day. 



TIME AND LENGTH SCALES 

The relevant nuclear burning time scale for the H-ingestion problem is the time scale for a proton to be captured by a 12 C: 

n= 12 

T ' ! P X( 12 C)pN a <av>„ C (p, 7 ) ' 



For the quantitative evaluation of the relevant time scales we use the pre-ingestion model at time fo shown in Fig. [2] (Sect. 3.2 1. The 
mass fraction of 12 C in that model is X( l2 C) = 0.36 and the density increases from p top = 1 .26 x 10 2 g/cm 3 at the top of the convection 
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FIG. B 13. — Time scales as a function of the mass coordinate in the convection zone, for proton capture by 12 C (blue solid line), as well as the MLT mixing time 
scale (green dash-dot) and the rate of reaction mixing time scale (red dashed) (see text for details). For this figure the tabulated reaction rate from (Angulo et al. 
[1999) was used. 



zone to pbot = 1-0410 x 10 4 g/cm 3 at the bottom of the convection zone. The nuclear reaction rate < av > 



2 C(p, 7 ) 



depends sensitively 



on the temperature which increases from T top = 2.2 x 10 7 K at the top to T\, ot = 2.9 x 10 8 K at the bottom of the convection zone. 



< av >a 



C(p,7) 



increases by 12 orders of magnitude accross the convection zone. 



The location of the peak H-burning due to H-ingestion takes place where the mixing time scale is the same as Ti2 C (p) (Ch. 4 Arnett 
|1996| l. The diffusion coefficient Dmlt for convective mixing is derived from the mixing-length theory (MLT). With an appropriate 
length scale / a mixing time scale can be obtained. For some properties the MLT mixing-length /mlt should be used: Zmlt = ^mljH p 
with aMLT = 1.7 the mixing-length paramter and H p the pressure scale height. This MLT mixing time scale is then tmlt = 'mlt/^mlt- 
As can be seen in Fig. jB 13 Ti2 C (p) = t mlt at m x = O.6O24M , a significantly larger mass coordinate than the location of the peak 
H-burning (~ 0.6005 M Q ) calculated in the stellar evolution model, as evident from the H-profile at t\ in Fig. [2] 

'mlt should not be used to estimate a mixing time scale relevant for rapid nuclear burning, since the rate of p-captures depends 
only indirectly on P. In fact, in the vicinity of the H-peak luminosity the pressure scale height is Hp ~ 1.4Mm which implies 
Imlt ~ 2.4Mm. This is much larger than the distance over which the rate of p-capture by 12 C (rate of reaction) increases significantly. 
It is this rate of reaction length scale that defines the width and location of the combustion flame for a given diffusion coefficient. A 
generalized length for any quantity <fi = (j>(r) may be defined as (Chapm an|1961 1 

1 

Ha, 



■ 



where is the rate of reaction length scale if we define (j> = pN a < av >i2 C ( p 7 ). The rate of reaction mixing time scale is then 
ta, = HI/Dmlt- As shown in Fig.B13 the mass coordinate where ri2 C (p) = ta, coincides very well with the location of peak H- 



burning (where as a result the mixing split occurs) at t = t \ in Fig. [2] 

At this location (m, ~ 0.6005 M©) the reaction length scale is Ha, ~ 330km which is the geometric scale of the flame that hydro- 
dynamic simulations including nuclear burn have to resolve. A simulation box that fits the An geometry of the entire He-shell flash 
convection zone needs to have a side length of 50Mm which corresponds to ~ 166 flame widths. In order to resolve the flame with 
at least 10 radial zones an aequidistant grid for a H-ingestion flash simulation needs to have a 1660 3 grid. 



